Dynamics of Si02-glasses 
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Results on the dynamics of silica are presented: vibrations and relaxations. Using molecular 
dynamics, glass structures are generated by rapidly quenching melts below the glass transition. For 
the local minima of the structures the vibrational density of states is determined, the structures of 
the eigenmodes are analyzed, and the influence of the single components is discussed. Relaxations 
are studied in both amorphous Si02 and silica melts. Our main focus is on the type of relaxation, 
i.e., whether the main contributions are caused by small atomic displacements, or by bond- breaking 
processes, e.g. creation and annihilation of dangling bonds. The eigenvector of the relaxations is 
found to be similar to the low-frequency one. 



PACS number(s): 61.43.Fs, 66.30. Fq 
I. INTRODUCTION 

Glasses belong to the oldest manufactured materials 
made by man. But despite the long traditions and abili- 
ties in producing glasses still sufficient information about 
the atomic structures in amorphous and glassy states is 
lacking. Due to the absence of periodicity found in crys- 
tals both experiment and theory are challenged to work 
with sophisticated methods. The technological relevance 
of ceramics has led to increased efforts to resolve the 
structures and to gain insight into the so-called structure- 
property-relationship which can be exploited for the de- 
sign of new materials. Neutron- or x-ray-scattering pro- 
vides information about the global arrangement of the 
atoms via the radial distribution of the pair-correlatiorffl. 
Electron spin resonance and positron annihilation is used 
to detect dangling bonds in the structures! From nu- 
clear magnetic resonance studies .one can get insight into 
the distribution of bond anglesa, however models are 
necessary to interpret the results. Due to the enor- 
mous increase in computer power, simulation has be- 
come a widely used tool to model amorphous structures. 
The applied techniques enco mpas s molecular dynamics 
(MDJTB, Mf4e Carlo (MCpliil and molecular mod- 
elling (MM)tLata, each of which has advantages and draw 
backs. Molecular modelling is successfully appliedr-tp co- 
valently bonded materials, e.g. SiO^tZl and GeOj^l. Its 
main advantage is that, without applying refined interac- 
tion potentials, high-quality structures can be modelled 
which show all the experimentally known data per con- 
structionem. The draw-back is the need of information 
about typical bond- lengths and bond-angles of structural 
units in the glass. Monte Carlo methods are applied 
to investigate the phase space and its properties-.and to 
study relaxations from a global point of viewta The 
Reverse Monte Carlo (RMC) technique is often used to 
build structures according ^experimental information, 
e.g. neutron-scattering dataEJ. The disadvantage of MC 
is that moveclasses may often be unphysical and do not 
give reliable insight into the microscopic dynamics of the 
systems. Molecular dynamics is widely used to construct 



models of amorphous state by rapidly quenching melts 
and analyzing the dynamics of the modelsEil. Due to the 
fact that high-quality interaction potentials are necessary 
to simulate both structure and dynamics of the materials 
under consideration the results of MD strongly depend 
on these potentials. In all cases (MM, MC and MD) the 
typical system sizes are of the order of several thousand 
atoms, which correspond to a volume of several ten thou- 
sands cubic Angstroms. These small system sizes cannot 
be used to simulate large scale effects, e.g. grain bound- 
aries (where a grain would be of the order of a typically 
available model size). 

Compared to crystalline phases amorphous materials 
show a richer dynamics which plays an important-^] 
with respect to low temperature thermal propertiescSI 
Apart from the long wavelength modes (sound waves) 
also present in crystals, in glasses strongly overdamped 
modes are observed at low frequencies. These strongly 
overdamped modes lead to a strong reduction of the ther- 
mal conductivity in glasses compared to their crystalline 
analogues. These modes can also be seen in the excess 
specific— heat at low temperatures . typically below a 
few K.EJ In the corresponding crystals one finds an in- 
creasepOijihe specific heat with T 3 according to the Debye 
modelEj'Ej. The low temperature behaviour of glasses is 
well described by the tunneling model assuming a con- 
stant distribution of two-level tunneling statesEZrEa. 

At somewhat higher temperatures th^ specific heat 
shows a b«jnp in c p /T 3 . Using RamarO and neutron 
scatteringtlJ one has deduced that thermally activated 
jumps over energy barriers and vibrations, which be- 
come anharmonic at the lowest frequencies, are responsi- 
ble for such a behaviour. The influence of these effects is 
well described by the soft potential model (SPM)EirE3 in 
which one assumes that the dynamics of soft regions in 
the structures is governed by soft potentials. The SPM 
is an extension of the standard tunneling model which 
is included in the SPM as its low temperature limit. 
The parameters describing the softness of both the two- 
and one- well potentials are distributed according to some 
probability distributions. Fitting these parameters to ex- 
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perimental data one finds that the regions comprising the 
active|-«jrbrating atoms are confined to about 20 to 100 
atomsHcJ. 

Recently, computer simulations of soft sphere gl asjas , 
(SSG) have supported the view of soft local modesc§E£l 
with entities comprising about 20 to 50 atoms. Simi- 
lar modes have since been observed iiu simulations of 



covalently bonded glasses, e.g.,— ■ SiC>2 nd SeE_l, in 



■.SicJ, and amor- 
3 and other sys- 



metallic glasses, such as Ni-ZrEj, P 
phous and quasi-crystalline Al-Za-M; 
terns, for example amorphous iceca. In an earlier simula- 
tion of amorphous silicon low-frequency localized—vibra- 
tions have been observed at coordination defectsc 3 ! 

Further and in addition to these periodic excitations, 
relaxations (aperiodic rearrangements) are observed in 
ultrasonic and dielectric relaxation experiments.C-3 These 
processes occur already at low temperatures [T > 5 K) 
and can be either tunneling processes or thermally acti- 
vated jumps over energy barriers. In the SPM local relax- 
ations are strongly correlated with local soft vibrations, 
assuming a smooth transition from vibrations to relax- 
ations. Experimentally this is supported by the similar- 
ity of the structure factors of both types of excitations! 3 !! 
From diffusion measurements in a metallic glass effective 
masseS|-of ten atoms have been derived for this collective 
motioned. Strong correlations between soft vibrations 
and hopping have been observed in moleculax. dynamics 
simulation of the SSGo and in amorphous SeeZl. Both re- 
versible and irreversible relaxations were observed which 
consist of a collective hopping of groups of atoms. Collec- 
tive hopping motions have also been observed in simula- 
tions of binary Lennard- Jones and soft-sphere mixtw-ps 
above and below the glass transition temperature^ JiS 
and in amorphous— Ar after introduction of vacanciesE9. 
Heuer and SilbeyEil systematically searched for double 
well potentials in a binary model glass at T = 0. They 
found a distribution of two well potentials in qualitative 
agreement with the soft potential model. 

In this paper we investigate the structure, vibrations 
and relaxations in amorphous silica. The following chap- 
ter is devoted to computational details. The analysis of 
the structures is given in section [II . The vibrations and 
related quantities are shown in section IV. The results 
of a detailed analysis of the ..modes similar to the one 
given by Taraskin and ElliottH are presented in section 
M. The study of relaxations in silica glasses and melts 
(section VI ) is the most important part of the work pre- 
sented here. To our knowledge this is the first exten- 
sive such study. As done in previous work on SSG and 
Se, the relaxations are analyzed with respect to jump- 
length, localization and element specific contributions to 
the atomic rearrangements. In the last section the results 
are discussed and conclusions are drawn. 



II. COMPUTATIONAL DETAILS 



A. Si02 and its potential 

Si02 exists in many different crystalline allotropes (e.g. 
a- and /3-quartz, high- and losfccristobalite, tridymite, 
keatite, coesite andpStishovite)!^. It is known to be a 
strong glass formerEltil. 

The atomic interaction potential used in our simula- 
tion was fitted by Vashishta et al.Q in order to reproduce 
structural and dynamical properties of both crystalline 
and glassy phases. The potential consists of two-body 
and three-body-terms. The two-particle interaction con- 
tains a long-range_Coulomb part which we treat by an 
Ewald summatiorJEj, and includes terms for steric repul- 
sion of the particles and polarizability of the atoms. 

The three-body-interaction is only relevant for Si-O- 
Si and O-Si-0 units, for the other possible combinations 
( Si-O-O, O-Si-Si, O-O-O, and Si-Si-Si) it is neglected. 
The three-particle terms favour the tetrahedral angle at 
the silicon (O-Si-O) and an angle Oo ~ 140° at the corner 
sharing oxygen (Si-O-Si). For details see Ref. @. 



B. Structures and Relaxations 

To produce silica glasses we quenched equilibrated 
melts below the glass transition temperature T g , the 
number density is kept constant during all MD-runs and 
corresponds to a density p — 2.20 g/cm 3 . An estimate 
of the glass-transition temperature was obtained by cal- 
culating the diffusion constant, D, of the system as a 
function of temperature using the relation 

D a = lim 1 <| R Q (0) - Ra(t) | 2 > , (1) 
t^oo bt 

where R Q (t) is the time-dependent position vector of a 
particle of type a € {Si, 0} and < ... > denotes the 
configurational average. 

In Fig. |l] the temperature dependent self diffusion con- 
stant for both silicon and oxygen shows a rapid drop at 
a temperature To rs 2500 K. This temperature To can 
be taken as a upper limit for the glass transition tem- 
perature T g , since below To the diffusive motion of the 
particles is effectively frozen out. In the high tempera- 
ture limit the ratio between the oxygen and silicon diffu- 
sion constant is about 1.2 which is less than one would 
expect from the mass ratio mgi/mo = 1-75 if the dif- 
ferent species diffuse according to their respective mo- 
mentum distribution. But this effect might be explained 
by the different "self" -interaction term in the Si-Si- and 
O-O-interaction. In comparison to experiment (where 
one measures T g w 1500 Ko) the calculated glass tran- 
sition temperature is much too high. This problem is 
common to all computer simulations and caused by too 
high quench rates, which exoeed the experimental ones 
by many orders of magnitudelia. 

We used constant quench rates of T < 3 • 10 13 K/s and 
quenched the glasses in two steps: first we quenched the 
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systems from T = 3500 K to 1700 K, monitored relax- 
ations and jumps in a temperature regime well above and 
below T g , and stored the atomic positions of the detected 
new minima. In a next step the temperature of these 
glasses are quenched with the same rate to 270 K. Again 
we monitored the relaxations and stored the atomic co- 
ordinates. We have generated 14 glasses with N = 1944 
atoms each. To study system size effects we have ad- 
ditionally constructed 10 smaller glasses with N = 576. 
The time step in our simulations is about 1.2 fs. Surface 
effects are reduced by applying periodic boundary condi- 
tions. During equilibration we used the NVT ensemble. 

The resulting structures are used to analyse relax- 
ations both around the glass transition temperature T g « 
2500 K and well below. To detect the relaxations in the 
course of the MD runs we monitored the atomic displace- 
ments, 



(2) 



of the potential ener 
ous unstable modest 
matrix are given by 
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prevents the occurrence of spuri- 
I. The elements of the dynamic 
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which are the mass-weighted second derivatives of the 
potential energy U with respect to the atomic positions. 
The frequency spectra are calculated from the frequencies 
of the 3iV — 3 vibrational modes a as 



Z(u) = 



1 



3^-3 



(4) 



where S is the discretized ^-function and (. . .) stands for 
the averaging over configurations. Due the limited sys- 
tem size the spectrum is cut-off at small frequencies. To 
close this gap we calculated the Debye spectrum 



where R n (t) is the position vector of atom n and R,\ 
gives its position in minimum i of the potential energy 
surface. If the total displacement of the atoms exceeded 
a cut-off value, and the residence time of the atoms in 
the new positions also exceeded a minimal period of at 
least three times the period of a typical soft vibrational 
mode, the new positions of the particles were accepted as 
a new minimum configuration. The cut-offs of displace- 
ment and resident time, respectively, are chosen to avoid 
spurious minima. Relaxations were observed by moni- 
toring the total displacement with respect to the initial 
minimum, see Eq. |2| the corresponding minima on the 
potential energy surface, which are identified by a final 
quench to T = K, were stored and analysed. 

To study the relaxations in glasses at elevated tem- 
peratures, the quenched glasses were heated to temper- 
atures between 270 K and 1670 K. At each temperature 
the glasses were observed for up to 250000 time-steps, 
corresponding to about 0.3 ns. To keep the tempera- 
ture constant we averaged T(ti) over a period of 20 time 
steps [T av = ^ Yli^m T(ti)]&nd scaled the velocities 
after each period by y/To/T av with To being the "de- 
sired" temperature and T av = the averaged temperature. 
This has the effect that potentially relaxing atoms are 
not "slowed down" by reducing their velocities. 



C. Vibrations 

After a final quench of the low temperature configu- 
rations to K by a steepest descent/conjugate gradi- 
ent algorithmic the minima of the potential hypersurfacc 
are identified. In these minima, we calculated the dy- 
namic matrix, diagonalized it and determined its eigen- 
values, the latter_correspond to the squares (ui 2 ) of the 
eigenfrequencies.EZI The numerically exact minimization 



^Dcbyc — 3 ^ 
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with 
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3N \ 1/3 
AwVJ 



(5) 



(6) 



and the average sound velocity c given in terms of the 
longitudinal and transverse velocities q and Ct- These 
arc calculated from the clastic constants of the glass by 
the usual relation q = \fcujp and c t = y c 4i/p where 
we employ the elastic isotropy of the glass. 

We calculated the elastic constants from the change in 
potential energy, AE under an applied strain 
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(7) 



1 e 1 p. 



a/3 



apjS 



a/37 



(8) 



Here the first term accounts for the work done against the 
forces for an ensemble which is not in equilibrium against 
volume changes where P a p is the virial of the forces. The 
third term is a correction for the volume change under a 
finite shear in such a lattice and the C a p 1 i are the elastic 
constants (en = Cim, c 44 = C 2 323)- 

We found for the glass the sound velocities q = 
6100m/s and c t = 4600m/s which are about 10% and 
20% higher than the experimental values. Taking the fi- 
nite size of the ensemble into account this implies that the 
lowest frequency sound wave which could be observed for 
the largest sample (N=1944 atoms) is about v w 1.5 THz. 
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III. STRUCTURE OF SILICA-GLASSES 

In Table [| the distribution of the coordination numbers 
CN for the single components at different temperatures is 
given. If the atomic distances m < 2.3 A (this distance 
corresponds to the minimum in the pair-correlation func- 
tion, see below) then atoms i and j are neighbours. One 
can clearly observe a reduction of the fraction of "coor- 
diantion defects", i.e. atoms with another coordination 
number than the most common one, with decreasing tem- 
perature. During cooling the local topology of the atoms 
improves due to relaxation of the structures. In the liq- 
uid phase at T = 3500 K nearly 10% of the particles are 
under- or over-coordinated. Whereas equal fractions of 
oxygen atoms are under- and over-coordinated, nearly no 
over-coordinated silicon is found. However, in the glassy 
state at T = 870 K the fractions of under-coordinated 
Si and O are reduced to less than 5% and less than 
3%, respectively. Over-coordinated oxygen atoms have 
nearly disappeared, and nearly no over-coordinated sili- 
con is found. In the glasses quenched from T = 270 K to 
T = K the total fraction of defects is even less than 3%. 
From experiments it is known that in silica coordination 
defects exist in which the silicon is only surrounded by 
three oxygen (oxygen vacancy )Ej the number of these de- 
fects is enhanced by irradiation. Our numbers of coordi- 
nation defects are comparable to other models and simu- 
lations. In molecular dynamics simulation of silica glasses, 
the numbers of defects can vary between less than 2%&E£l 
and 6% to 8% as pointed out in Ref. ^| and references 
therein, the effective coordination numbecs strongly de- 
pend-on the applied interaction potentialsa and quench- 
ratesEj. Using molecular modelling techniques vitreous 
structures without any defects can be constructed per 
definitioner r&J&B Typically these structures are cluster 
models, i.e. without periodic boundary conditions. 

Quantities well-suited to gain insight into structural 
correlations and comparing simulations and experiments 
are the total and partial pair-correlation functions. We 
take the usual definition of the partial pair-distribution 
functions g a f}(r) 



< n a p(r) > Ar = Airr 2 Arp N cpg a i3{r) 



(9) 



where n„^(r)Ar is the number of particles of species /3 in 
a shell of thickness Ar and radius r around a particle of 
species a and < ... > denotes the configurational average. 
Pm = N/V is the number density and cp = Np/N the 
concentration of species (3. The partial pair-correlation 
functions obtained at T = 270 K shown in Fig. ||(b- 
d) are used to assign the peaks in the total g(r), Fig. 
§(a), to the partial contributions. From the integral over 
the first maximum of the pair-correlation functions one 
can deduce average coordination numbers: each silicon 
has 3.96 oxygen neighbors. The coordination numbers 
of the O-O-correlation reaches a number 6 for the range 
from 0.25 nm to 0.3 nm. For Si-Si we observe a value 
of 4 in the range from 0.28 nm to 0.35 nm. From these 



partial pair-correlations one can deduce that the peak at 
0.42 nm with a shoulder at about 0.35 nm is a superposi- 
tion of the second peak in the Si-O-contribution and the 
shoulder of the second O-O-peak. The fifth peak in the 
pair-distribution function at 0.5 nm in Fig. | appears to 
be a superposition of the second peaks of the O-O- and 
Si-Si-pair correlations. 

In Fig. | the pair-distribution function d(r) = 
AirrpN{g{r) — 1) averaged over an ensemble of 14 glasses 
with N = 1944 atoms is shown. Togethep,with the re- 
sult of an electron scattering experimented. Up to For 
r s» 0.6 nm, i.e. twice the average Si-Si-distance, both 
curves agree very well. (The wiggles below the first peak 
in the experimental curve are an artefact.) 

The bond-angle distributions are plotted in Fig. |J. 
The O-Si-0 distribution has its peak-position at 109° 
and a full width at half maximum (FWHM) of 10°. Such 
a narrow distribution of the tetrahedron-angle suggests 
slightly distorted tetrahedra as basic structural unit. The 
Si-O-Si distribution shows a peak at 145° and a FWHM 
of 25°. This broader angle distribution refelcts the disor- 
der in the silica glasses. In the O-O-Si angle distribution 
the peak at 35° stems from atoms belonging to the same 
tetrahedron. It is directly connected to the sharp peak in 
the O-Si-0 angle distribution. The broad distribution be- 
tween 80° and 180° stems from configurations where the 
two O-atoms have different Si's as nearest neighbours. In 
the Si-Si-0 angle distribution the peak around 20° stems 
from two silicon atoms bonded to the same oxygen. This 
peak is directly connected to the Si-O-Si distribution. If 
an oxygen is not a bridging atom between two silicon 
atoms, but instead neighbouring the adatom, we find an 
angle between 60° and 160°. The strong peak at 60° 
in the 0-0-0 angle distribution is from triplets of O's 
bonded to the same Si-atom. For a situation in which 
the oxygen atoms are bonded to different Si the 0-0-0 
angles are always larger. For the Si-Si-Si angle we find 
a broad distribution between 80° and 180° with a maxi- 
mum at about 105°. This might give a hint to patterns in 
which the silicons are surrounded by strongly distorted 
"Si-tetrahedra" . The distance between the next-nearest 
silicons in the triple is about 0.5 nm. In triples where 
all Si are at their nearest neighbor distance (w 0.3 nm) 
we find Si-Si-Si angles of about 60° which explains the 
small vicinal maximum in the distribution. This points 
to the existence of 6-membered Si-O-rings in silica. In 
all, the angle distribution clearly reveals corner-sharing 
tetrahedra of Si04. 



IV. SPECTRA AND LOCALIZATION OF 
VIBRATIONAL MODES 

Using Eq. [| we show the vibrational density of states 
as Z(y) in Fig. |^ averaged over all configurations to- 
gether with the Debye spectrum, Eq. ||. One sees a clear 
enhancement of the glassy spectrum at low frequencies 
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which reaches well beyond frequencies where system size 
could be of importance. Reflecting some shortcomings 
of the interaction potential thepfixperimentally observed 
double peak at high frcqucnciesEj cannot be reproduced. 
The Bose-peak, i.e. the maximum of Z{y)/v 2 , lies at ~ 
1.6 THz high©* than experiment where one observes it at 
about 1 THzEa. From the density of states we can calcu- 
late the vibrational specific heat at constant volume c v . 
In the harmonic approximation one has per atom 



c v — 3k b / du> 



h \ 2 

/sinh 2 (^/2fc B T) 

2Kb ) 



(10) 



The vibrational specific heat of a perfect crystal at low 
temperatures is oc T 3 . It is, therefore, usual practice to 
plot c v /T 3 . In such a plot a Debye spectrum gives a con- 
stant whereas in glasses an increase above this constant 
is found. 

Fig. ^| shows this behavior for the silica glasses. The 
solid line shows the values gained from the spectrum 
of Fig. 1 which is corrected for the finite size of the 
simulated glass by adding a Debye contribution up to 
a frequency smaller than the lowest possible phonon 
frequency. This correction amounts to a fraction of 
7.7 x 10~ 3 of all modes. The resulting values for the 
specific heat are shown by the full-line. The diamonds 
corresponds to experimental valuesEa. The too high val- 
ues of the sound velocities in our model cause a too high 
value of vd which leads to a Debye contribution to cy 
which is too small by a factor of f=s 1.6. 

The phonon eigenvectors describe the structure of the 
vibration and can, e.g. be used to determine the degree 
of the localization of the vibration. There are two com- 
monly used measures of localization, the effective mass 
and the participation raticOE3. The effective mass is 
given in terms of the eigenvector as 



TOeff(cr) = m/le 1 ^)! 5 



(11) 



Here we have assumed that the 3./V-dimensional unit vec- 
tor of mode a is normalized and e n (a) stands for the 
vector formed from the three components on atom n. 
Atom number 1 is chosen as the atom with the largest 
displacement. m e s/m is a measure for the number of 
atoms which effectively carry the kinetic energy of the 
vibrational mode. This definition is limited to small sys- 
tem sizes when the long range tails of the modes are not 
too important. In the following we will mainly use the 
participation ratio: 



N 

p{a)= ( iVX|e»| 



(12) 



For a translation one has p — 1 and for a vibration of a 
single atom with all others at rest p = 1/N. This scal- 
ing with 1/N should hold for all localized modes. Fig. 



fj] shows the participation ratios for the two system sizes 
with N — 576 and N — 1944 atoms. The high-frequency 
modes are highly localized and their participation ratios 
show the expected scaling with system size. At the low- 
frequency end we observe (quasi-)localized modes. The 
localization corresponds to clusters comprising about 10 
to 50 atoms. The strongest localized modes have partici- 
pation ratios of 0.012 corresponding to an entity compris- 
ing about 8 atomSj—These findings are in agreement with 
ones by Jin et alEj and Taraskin and Elliotto. How- 
ever, their participation ratios do not show the expected 
1/N scaling. This will be due to interaction effects be- 
tween the modes in the larger structures which results 
in an increase of the participation ratio. In the larger 
systems there are more of similar frequency, and in par- 
ticular more low frequency extended phonons which will 
mix with the quasi-localized modes and raise their par- 
ticipation ratio. Additionally this causes an increased 
interaction between quasi localized modes of similar en- 
ergy, since due to de-localization these modes lose their 
"splendid isolation" and "feel" also other localized vi- 
brations, too. Due to this interaction of quasi localized 
modes their participation ratio is increased. For a finite 
concentration c of interacting modes the scaling factor 
1 /N in Eq. [l2] will be replaced bx c. Such effects have 
also been found in a model glasst 3 ] where the modes are 
deconvoluted such as to reconstruct the "naked" contri- 
butions: phonons and localized vibrations. 



V. ANALYSIS OF VIBRATIONAL MODES 

To gain more insight into the dynamics of the vibra- 
tional modes we analysed the motion of the single com- 
ponents and their contribution to the eigenvectors. The 
averaged contribution of the oxygen to the mass- weighted 
eigenmodes in comparison with the silicon is depicted in 
Fig. |8|. This contribution is given as the mean square 
displacement of each component a G {Si,0} to the eigen- 
vectors 



Z a (u) 



1 



3N- 3 ^ 



i=l 



(13) 



with N a the number of the Si- and O-atoms, respectively. 
Just below the high-frequency peak the silicon atoms 
contribute stronger to the eigenvectors than the oxygen 
atoms, whereas the eigenmodes of the low-frequency part 
of the spectrum and the high-frequency modes are clearly 
dominated by the oxygen motion. To learn about typical 
motions we calculated the angle between the displace- 
ment of atom i in the modes (ej), and the bonds of this 
atom to its nearest neighbors j, i.e. rbonds = ry 



bonds 



arccos 



^i^bonds 

led I rbonds I 



(14) 



In Fig. H we show the weighted distribution of the 
angles for the low-frequency modes of each glass with 
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respect to the components Si and O. We weighted with 
e 2 in order to suppress marginal contributions of only 
slightly moving atoms. The contributions of both ele- 
ments are peaked at 90° which means the motion of both 
Si and O is mainly perpendicular to the bonds. In the 
case of oxygen this peak is very sharp, whereas the sili- 
con atoms have stronger contributions to either smaller 
or larger angles, this indicates that the silicon atoms have 
contributions parallel tOj_their bonds. Following the work 
by Taraskin and ElliottEj we calculated the projections 
of the vibrations of the structural subunits, such as the 
Si04-tetrahedra and the SiOSi-units, onto the symmetry 
vibrational coordinates, i.e. internal vibrations, such as 
symmetric and asymmetric stretching, bending modes, 
and rigid body motions, namely rotations and trans- 
lations. The translations of the sub-units are neglible 
(< 1%) in comparison to rotations and internal vibra- 
tions. To calculate the internal motions the center of 
these subunits are fixed, i.e. the vibrations of the sur- 
rounding atoms are shifted relative to the central atoms. 
The remaining motions can be either rigid rotations or 
internal vibrations, such as stretching or bending modes. 
In Fig. [l^ we show the contributions of the projections 
onto the SiOSi subunits to the total DOS. The modes in 
the low-frequency range have a strong rotational charac- 
ter, whereas in the high-frequency peak an asymmetric 
stretching mode describes well the typical vibrations. In 
the mid part of the spectrum symmetric bending and 
with less extension also symmetric stretching contribute 
to the modes. Projecting the modes onto the symme- 
try vibrations of the SiC>4 structural subunits (see Fig. 
|ll| ) we find the high-frequency peak to be described well 
by both asymmetric (strong contribution) and symmet- 
ric (weak contribution) stretching modes. Sipee in the 
experimental DOS one observes a double peakEJ, we can 
speculate that our interaction potential gives a too low 
concentration of symmetric stretching modes and a too 
high concentration of asymmetric modes. The center of 
the spectrum has contributions from bending modes with 
F2 and E symmetry. The low-frequency side of the spec- 
trum is again dominated by rotations. This is in agree- 
ment with the observation that the atoms move perpen- 
dicularly to their bonds. 



VI. RELAXATIONS ABOVE AND BELOW THE 
GLASS TRANSITION TEMPERATURE 

At the start of the simulations the temperature of the 
configurations is about 3500 K. Cooling with a quench 
rate of about T < 3 • 10 13 — the structures cooled to 
T = 1700 K after 50000 time steps w 60 ps. During cool- 
ing we observe initially a rapid succession of strong relax- 
ations. After about 20 ps and at a temperature of circa 
2800 K ( slightly above T g ), the activity of the structures 
slows down, i.e. jumps become less frequent and the sys- 
tems become temporarily trapped in metastable regions 
of configuration space. 



During cooling several new local minima of the poten- 
tial energy surface are visited by the "MD-walker" via 
jumps over energy barriers. To determine the localiza- 
tion of the hopping processes and to test the SPM as- 
sumptions, we calculated similar to Eq. [n] the effective 
mass of the relaxation: 



(AR) 2 
IARLJ 



(15) 



where (AR) 2 = 5^ n (RJ, — ^-n) 2 * s the square of the total 
distance AR between two successive minimum configura- 
tions (called "initial" and "final" positions of the jump). 
R^* denotes the respective initial and final positions of 
atom n, jAR^J is the maximal distance a single atom 
jumps in this relaxation, and m max is the mass of the far- 
thest jumping atom. Similarly to Eq. [l^ we calculated 
the participation ratio par to determine the localization 
of the hopping processes, 



PAR 



AR 4 



(16) 



where R^ and R£ denote the initial and final position 
of atom n in the relaxation and AR is the total jump- 
length. The participation ratio has the value n/N if n 
atoms are involved in the hopping process, and par. = 1 
if all atoms contribute equally to the relaxation. We find 
that the participation ratio of jumpsgrows roughly lin- 
early with the jump distance, Fig. O. Comparing the 
results for the two temperature intervals one clearly ob- 
serves an increase of both AR and par with tempera- 
ture. The gap between the two temperature regimes is 
an artifact of the cut-off procedure used for the higher 
temperatures. 

To study relaxations in a "metastable equilibrium" we 
heated the quenched glasses to T = 270, 870 and 1670 K. 
For the jumps of the glasses which are heated from T = 
K to the desired temperature we monitored again the 
minima visited in course of the MD runs and determined 
the corresponding jumps and their participation ratios. 

As example we plotted the time evolution of the energy 
and the displacement in Fig. [l4| for the glasses at T = 
870 K. In the low-temperature regime we can follow an 
overall relaxation of the structures towards more stable 
configurations. This trend becomes even stronger in the 
higher temperature regime. The ensemble-averaged par- 
ticipation ratios for the equilibrated temperatures (equili- 
brated during heating) show a slight temperature depen- 
dence, par = 0.020 ± 0.008, 0.057 ± 0.035, 0.061 ± 0.036 
for T = 270, 870, and 1670 K, respectively. Quite in 
contrast though, the jump-lengths increase linearly with 
temperature, and we find AR = 1.81 ± 0.75, 4.72 ± 1.38, 
and 8.27 ± 2.49A for T = 270, 870, and 1670 K, respec- 
tively. The typical "resident time" per minimum is about 
15 ps. In studies of amorphous Se the resident times per 
minimum is found to be about 30 ps at 0.05 T g O 

To identify active regions we calculated a correlation 
between the observed relaxations 
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CRRi 



J2 n AR n AR' n 
ARAR' 



(17) 



Since we calculated the product of the absolute values 
AR and not the scalar product between the two jump- 
vectors AR and AR , we also observed active regions 
where almost the same atoms contributing to a relax- 
ation AR jump perpendicularly to the previous relax- 
ation AR. In that case crri is close to 1. If the jumps 
are uncorrelated crri becomes n/N , with n the typical 
number of atoms contributing to the relaxations, and N 
the total number of atoms in the configuration. 

Calculation of the correlation between successive 
jumps (see Eq. 0) at the lowest considered temperature 
(T = 270 K) yields values of 0.585±0.314. Astonish- 
ingly, for the glasses at T — 870 and 1670 K very similar- 
values (0.654±0.089 and 0.666±0.082, respectively) are 
observed. The broader distribution of values in the low- 
temperature simulation may be explained by the obser- 
vation that at the lowest temperature successive relax- 
ations are either reversible jumps (high correlations) or 
jumps where different parts of the active regions of the 
systems contribute to the relaxations at different times 
(low correlations). At higher T these local regions aggre- 
gate to a larger complex where the single entities are not 
resolved in the jump-process. This can explain the rel- 
atively high correlations between the jumps at higher T 
where several "active" regions in the glasses jump at the 
same time. However, at the higher temperatures succes- 
sive jumps are not reversible, but the same (local) regions 
are activated throughout the observation time. 

To see whether there exists a typical "relaxation- 
mode" we calculated the projection of the relaxations 
observed at T = 870 K and 1670 K onto the bonds of the 
atoms, Fig. n5[ 



Oi ri bonds = arccos 



ARjTbonds 

IarTk 



onds 



(18) 



where AR^ is the relaxation vector of atom i. The 
distribution is Gaussian in the case of the oxygen, and 
centered around 90° indicating that the O-atoms relax 
mainly perpendicularly to the bonds. For the silicon 
atoms we observe a peak at low angles which reflects 
the fact that the Si-atoms relax mainly in direction of 
one bond. Since a silicon is surrounded tetrahedrally by 
the O-atoms, the projection onto the remaining bonds 
will give the tetrahedron angle at about 109°, which is in 
satisfying agreement with the peak at about 105°. Note 
that the ratio between the peaks at 5° and 105° is about 
1:3 which is the ratio one would expect in a tetrahedron 
if the center atom moves parallel to one of its bonds. 
In analogy to the mode-analysis we calculated the pro- 
jections of the relaxations onto the symmetry modes of 
SiOSi- and SiOvsubunits. Since in this study it has 
been important to consider the "real" contributions of 
the atoms and not only the "relative" displacements as 
in the case of the eigenmodes we do not subtract the 
displacement vector of the central atom. In the case of 



the SiOvsubunits and relaxations at T = 870 K and 
1670 K we do find nearly no translations of the sub- 
unit worth mentioning and negligible contributions for 
symmetric stretchings of the oxygen atoms; 0.21 of the 
relaxations of the oxygen atoms stem from asymmetric 
stretching, 0.40 of the jumps are due to bendings in F2- 
symmetry, 0.07 are bending with E-symmetry and 0.32 
of the oxygen motions come from rotations. Consider- 
ing the SiOSi-subunits we do not record a considerable 
amount of translations of the SiOSi-pattern; and only 
0.037 of the Si-relaxations are due to symmetric stretch- 
ing. The contribution of asymmetric stretching is about 
0.26, symmetric bendings contribute at about 0.23 to the 
jumps, and 0.47 come from rotations. The strong contri- 
butions from bending and rotations for Si and O are in 
agreement with the picture that the typical relaxation is 
perpendicular to the bonds. For the relaxations of sili- 
cons (in the SiOSi-subunits) we find that about 30% of 
the relaxations can be ascribed to stretching (3.7% sym- 
metric and 26% asymmetric) , for the oxygen relaxations 
in the SiOvunits only 21% of the relaxations can be de- 
scribed by asymmetric stretching. 

The investigation of relaxations with respect to bond- 
changes (creation and annihilation) , where a bond is cut 
off at 2.3 A, at constant temperatures shows that at 
T = 270 K the main contribution to a relaxation is due 
to small changes of the positions of the atoms, bond- 
breaking occurs only in a few cases. At T = 870 K about 
0.1% of the bond situation is altered by relaxations from 
one minimum to another. At the highest temperature 
the number of bonds change up to 0.3% per jump. 

During the cooling runs we stored configurations for 
some of the glasses at T = 870 K and T = 270 K. After 
an equilibration time of 60 ps we started relaxation runs 
for these "young glasses" at T = 870 K and T — 270 K, 
respectively. Again relaxations were monitored for about 
0.3 ns. 

The participation ratios averaged over the ensembles 
for the equilibrated temperatures T — 270 K (4 glasses) 
and T = 870 K (2 glasses) are p A R = 0.020 ± 0.011 and 
Pat?. = 0.047 ± 0.030, respectively. The jump-lengths 
increase with temperature, we find AR = 2.05 ± 0.77A, 
5.86±2.15A again for T — 270 K and 870 K, respectively. 
Now the typical "resident time" per minimum is about 
10 ps. Note that the jump-lengths of the "young" glasses 
are larger and the resident time is smaller than in the 
corresponding "older" glasses. These results although 
based on weaker statistics may point to the fact that the 
"young glasses" are more active in terms of longer jump- 
distances and shorter relaxation-times than those glasses 
aged for a longer time. 



VII. DISCUSSION AND CONCLUSION 

We have performed molecular dynamics calculations 
on silica glasses and melts. The amorphous structures 



7 



were generated by rapid quenches of melts. The analysis 
of the glasses show a satisfying agreement with exper- 
imental results for the pair-correlation functions. The 
angle-distributions reveal the existence of a network of 
corner-sharing tetrahedra. During cooling we observed 
an appreciable optimization of the local topology of the 
structures. At T = 3500 K the number of dangling bonds 
of Si and O is about 9.4% and 4.5%, respectively. Caused 
by structural relaxations the number of these defects de- 
creases to 3.9% and 2.2%, respectively, at T = K. This 
improvement can also be seen in the shift of the peak in 
the Si-Si pair-correlation function from less than 0.3 nm 
in the liquid, in which the Si-Si-distance appears as a 
shoulder of the O-O-peak, to 0.31 nm (at T = 270 K). 
The Si-Si-distance becomes visible as a separate peak in 
the total pair-correlation function. 

One main focus of our study concerns the dynamics of 
the structures with respect to both periodic motions (vi- 
brations) and aperiodic ones (relaxations). Calculation 
of the element specific contributions shows the silicon mo- 
tion to be dominant just below the high frequency modes 
(about 28 THz). As to the rest of the spectra, the oxygen 
vibrations become the main contribution to the mode. 
The participation ratios show that all modes in the high- 
frequency regime are localized, comprising entities with 
less than 15 atoms. At the lowest frequencies the modes 
— about 0.2% of the total spectrum — are typically quasi- 
localized vibrations with participation ratios of less than 
0.22 for N — 1944 atoms. The effective numbers of atoms 
contributing to the quasi-localized modes range from 10 
to 50 particles. The atoms strongly participating in those 
low-frequency modes move preferentially perpendicularly 
to their bonds. Investigating the subunits of the struc- 
tures such as Si04-tetrahedra and "SiOSi-fragments" we 
find that the motion of the low-frequency modes can be 
described by rotations (most important feature) with an 
addition of bending type modes. This picture is in agree- 
ment with the one of rotations of coupled tetrahedraE3. 
The high-frequency vibrations can be described well by 
stretching modes. 

We studied relaxations of Si02 structures in the melt 
(during cooling) and in the glass at several temperatures 
both during cooling and "equilibration" . During cooling 
from 3500 to 1700 K we observed rearrangements of the 
structures with participation ratios of less than 0.3 (for 
N = 1944 atoms) connected to effective masses ranging 
from 20 to 150 which are somewhere between clear local 
relaxations and global changes. This observation indi- 
cates that several regions of the structure are active at 
the same time. Lowering the temperature (from 1700 K 
to 300 K) the picture of locally changing structures be- 
comes clearer and the participation ratio of the jumps 
is less than 0.1. Analysis of the relaxations observed 
during the "equilibration" runs shows that the jumps at 
low temperature (T = 270 K) typically lead to small 
changes of the atomic positions with nearly no bond- 
breaking processes. At higher T we observe bond-changes 
of the atoms (both bond-breaking and bond-creations): 



At T = 870 K the numbers of bonds changed by about 
1 bond out of 1000, and at T = 1670 K of up to 3 bonds 
out of 1000 changes per relaxation. These relaxations are 
local processes comprising typically 10-50 atoms. For the 
lowest temperature we find about 10 atoms contributing 
strongly to the relaxations whereas at higher T (870 and 
1670 K) about 20-50 atoms are contributing. However, 
the jump-length increases clearly linearly with tempera- 
ture. In the course of the MD simulations the structures 
relax to energetically lower lying minima. During aging 
of the structures this trend seems to slow down. Further 
corroboration of this lies in the observation of longer res- 
idence times of the configurations in each minimum and 
smaller jump-lengths towards the end of our simulations. 

For all temperatures below the glass-transition we ob- 
served correlated relaxations, i.e. several active regions 
coupled via a few atoms contributing to the jumps. At 
the lowest temperature several active regions can be dis- 
tinguished whereas these regions are "coupled" at higher 
T, i.e. the atoma_pomprising these regions are active at 
the same tirndfj'tJ. 

Analyzing motion patterns of the relaxations has 
shown that the oxygen atoms contribute dominantly to 
the relaxations by typically small displacements mainly 
perpendicular to the bonds. Nearly half of the oxygen 
motion is due to bond bendings and 30% to rotations 
(0.32). The silicon atoms exhibit also jumps more closely 
in the direction of one of its oxygen neighbours. This 
shows a more stretching-like motion of the Si-relaxation. 
However, the total contribution of this motion pattern 
amounts about only 0.3. The rest of the silicon mo- 
tion can be attributed to rotations (0.47) and symmet- 
ric bendings (0.23). Comparing the motion patterns of 
the relaxations to the eigenmodes (vibrations) the obvi- 
ous similarity between the motion pattern found in both 
low-frequency vibrations and relaxations (main contribu- 
tions are from rotations and bending) indicates a strong 
connection between low-frequency modes and local re- 
laxations and supports the assumption of the soft po- 
tential model. Such a correlation between local relax- 
ations and-eigenmodes has also been found in amorphous 
sclcniumJ£3 However, the relaxations have stronger con- 
tributions from bending- and stretching-motions than the 
low-frequency modes which lead us to the conclusion that 
also higher-frequency modes are necessary to describe re- 
laxations and are acting as impetus for localized jump- 
processes. 
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TABLE I. Distributions of the coordination numbers 
(CN) for both silicon and oxygen at T = 3500 K, 870 K and 
K, above and well below the glass transition temperature, 
respectively. 



10 



FIG. 1. Element specific self diffusion constants in SiC>2 for 
Si (o) and O (+) plotted versus temperature. 

FIG. 2. Total pair correlation function g(r) (a) and the 
partial pair correlation functions (b-d) in amorphous SiC>2 at 
T = 270 K. 

FIG. 3. Comparison between the calculated pair distribu- 
tion function d(r) and the experimental-result by D. Heine- 
mann using electron scattering (see Ref.c9). 



FIG. 4. Bond angle distribution in silica glasses. The 
O-Si-0 distribution has the peak-position at 109° and a 
FWHM of 10° . The Si-O-Si distribution has the peak-position 
at 145° and a FWHM of 25°. 

FIG. 5. Density of states Z(v) vs. v obtained by calculation 
of the Fourier transform of the displacement autocorrelation 
function using the equation of motion (EOM) method (solid 
line) and diagonalization (dashed line). The dotted line dis- 
plays the Debye DOS. 

FIG. 6. Specific heat cv as cv /T 3 in units of [fiJ/gK 4 ] vs. 
T in K. Plot is double-logarithmic. Solid line is the specific 
heat with Debye-correction and dotted line without correc- 
tion. Diamonds are experimental values of specific heat c p by 
Zeller and Pohl. 

FIG. 7. Participation ratio of the vibrational eigenmodes 
for glasses with N = 1944 (o) and TV = 576 (+) atoms. 

FIG. 8. Averaged element specific contributions to the vi- 
brational spectrum for 14 glasses with N = 1944 atoms. 
Total DOS (solid line), O-contribution (dashed line) and 
Si-contribution (broken line) are plotted versus frequency in 
THz. 

FIG. 9. Weighted averaged element specific distribution of 
the angles a e bonds between the atomic eigenvector and the 
bonds for all low-frequency modes. The full line show the an- 
gle-distribution found for silicon, the broken line for oxygen. 

FIG. 10. The vibrational DOS (solid line) and the 
partial contributions for the projections of the relative 
atomic displacement onto the symmetry coordinates of 
the SiOSi-subunits with the contributions from rotations 
(dashed-dotted line), symmetric stretching (broken line), 
asymmetric stretching (dashed line) and from symmetric 
bending (dotted line) plotted versus frequency in THz. 



FIG. 11. The DOS (solid line) and the partial contribu- 
tions for the projections of the relative atomic displacements 
onto the symmetry coordinates of the Si04-subunits with the 
contributions from rotations (dashed-dotted line), symmet- 
ric stretching (broken line), asymmetric stretching (dashed 
line) and from bending in E-symmetry (dotted line) and 
in F2-symmetry (dashed-dot-dotted line) plotted versus fre- 
quency in THz. 

FIG. 12. Displacement of O vs. displacement of Si in single 
jumps. The line shows the calculated O displacement if one 
scales AR Sl by the square root of the mass ratio O to Si and 
by the ratio of the atomic numbers. 

FIG. 13. Participation ratio for the quenching regimes vs. 
displacement. The line is a guide to the eye. 

FIG. 14. Change of potential energies per atom averaged 
over an ensemble of 8 MD-runs at T = 870 K (a). Total dis- 
placements AR 2 over an ensemble of 8 MD-runs at T — 870 K 
(b). 



FIG. 15. Distribution of the angle a r b on ds between the re- 
laxations r and the bonds of the atoms, full line Si-relaxation, 
dashed line O-relaxation. 
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